阵列信号DOA估计系列(二).导向矢量与空间FFT(附代码)

在DOA估计里面，经常会看到导向矢量这个名词，也有的地方叫方向矢量，方向矩阵，基本上都是array steering vector 的翻译。

本文首先对均匀线阵（ULA，uniform linear array）的导向矢量做一个推导说明，然后介绍一个最基本的DOA估计方法——空间FFT。

在上一篇文章中（ 阵列信号DOA估计系列(一).概述.），已经对空间相位差做出了详细说明，下面在空间相位差的基础上，引出均匀线阵的导向矢量。

依然来看这样一个阵列：在一条直线上，均匀排布着 M M M 个阵元 假设远场信号为 s ( t ) s(t) s(t),阵列以第一个阵元（上图编号为0）为参考，则整个阵列接收到的信号为 x ( n ) = [ s ( t ) , s ( t ) e − j 2 π f 0 d s i n ( θ ) c , s ( t ) e − j 2 π f 0 2 d s i n ( θ ) c , s ( t ) e − j 2 π f 0 3 d s i n ( θ ) c , ⋯   , s ( t ) e − j 2 π f 0 ( M − 1 ) d s i n ( θ ) c ] \mathbf{x}(n) = [s(t),s(t)e^{-j2\pi f_0\frac{dsin(\theta)}{c}},s(t)e^{-j2\pi f_0\frac{2dsin(\theta)}{c}},s(t)e^{-j2\pi f_0\frac{3dsin(\theta)}{c}},\cdots,s(t)e^{-j2\pi f_0\frac{(M-1)dsin(\theta)}{c}}] x(n)=[s(t),s(t)e−j2πf0​cdsin(θ)​,s(t)e−j2πf0​c2dsin(θ)​,s(t)e−j2πf0​c3dsin(θ)​,⋯,s(t)e−j2πf0​c(M−1)dsin(θ)​]稍作变形，即可得到 x ( n ) = [ 1 , e − j 2 π f 0 d s i n ( θ ) c , e − j 2 π f 0 2 d s i n ( θ ) c , e − j 2 π f 0 3 d s i n ( θ ) c , ⋯   , e − j 2 π f 0 ( M − 1 ) d s i n ( θ ) c ] s ( n ) \mathbf{x}(n) = [1,e^{-j2\pi f_0\frac{dsin(\theta)}{c}},e^{-j2\pi f_0\frac{2dsin(\theta)}{c}},e^{-j2\pi f_0\frac{3dsin(\theta)}{c}},\cdots,e^{-j2\pi f_0\frac{(M-1)dsin(\theta)}{c}}]s(n) x(n)=[1,e−j2πf0​cdsin(θ)​,e−j2πf0​c2dsin(θ)​,e−j2πf0​c3dsin(θ)​,⋯,e−j2πf0​c(M−1)dsin(θ)​]s(n) 可以看到，接收到的信号向量 x ( n ) \mathbf{x}(n) x(n)是一个向量 [ 1 , e − j 2 π f 0 d s i n ( θ ) c , e − j 2 π f 0 2 d s i n ( θ ) c , ⋯   , e − j 2 π f 0 ( M − 1 ) d s i n ( θ ) c ] [1,e^{-j2\pi f_0\frac{dsin(\theta)}{c}},e^{-j2\pi f_0\frac{2dsin(\theta)}{c}},\cdots,e^{-j2\pi f_0\frac{(M-1)dsin(\theta)}{c}}] [1,e−j2πf0​cdsin(θ)​,e−j2πf0​c2dsin(θ)​,⋯,e−j2πf0​c(M−1)dsin(θ)​] 乘以一个标量 s ( n ) s(n) s(n)，且此向量是信号来波方向 θ \theta θ 的函数。

做如下定义： a ( θ ) ≜ [ 1 , e − j 2 π f 0 d s i n ( θ ) c , e − j 2 π f 0 2 d s i n ( θ ) c , ⋯   , e − j 2 π f 0 ( M − 1 ) d s i n ( θ ) c ] \mathbf{a(\theta)}\triangleq[1,e^{-j2\pi f_0\frac{dsin(\theta)}{c}},e^{-j2\pi f_0\frac{2dsin(\theta)}{c}},\cdots,e^{-j2\pi f_0\frac{(M-1)dsin(\theta)}{c}}] a(θ)≜[1,e−j2πf0​cdsin(θ)​,e−j2πf0​c2dsin(θ)​,⋯,e−j2πf0​c(M−1)dsin(θ)​] 于是，接收到的信号可以表示为 x ( n ) = a ( θ ) s ( n ) \mathbf{x}(n) = \mathbf{a(\theta)}s(n) x(n)=a(θ)s(n)可见，向量 a ( θ ) \mathbf{a(\theta)} a(θ) 将 阵列信号DOA估计系列(一).概述 所述空间相位查全部包含在内，即包含了信号 s ( n ) s(n) s(n) 的角度信息。这就是我们常说的导向矢量。

这里，可以看到以下几点信息：

为方便描述，后面的文章中，均以第一个阵元为参考，并把导向矢量和接收信号写为列向量。 同时，为帮助读者理解，本节的向量和矩阵均标注了维数

当有一个信号 s ( n ) s(n) s(n) 从 θ \theta θ 方向入射到阵列时，接收信号可以表述为 x ( n ) M × 1 = a ( θ ) M × 1 s ( n ) \mathbf{x}(n) _{_{M\times 1}}= \mathbf{a(\theta)_{_{M\times 1}}}s(n) x(n)M×1​​=a(θ)M×1​​s(n)

当有 N N N 个信号 s 1 ( n ) , s 2 ( n ) , ⋯   , s N ( n ) s_{_1}(n),s_{_2}(n),\cdots,s_{_N}(n) s1​​(n),s2​​(n),⋯,sN​​(n) 分别从 θ 1 , θ 2 , ⋯   , θ N \theta_1,\theta_2,\cdots,\theta_N θ1​,θ2​,⋯,θN​ 入射到阵列时，按照叠加的思维，收信号可以表述为 x ( n ) M × 1 = a ( θ 1 ) s 1 ( n ) + a ( θ 2 ) s 2 ( n ) + ⋯ + a ( θ N ) s N ( n ) \mathbf{x}(n)_{_{M\times 1}} = \mathbf{a(\theta_1)}s_{_1}(n)+ \mathbf{a(\theta_2)}s_{_2}(n)+\cdots+ \mathbf{a(\theta_N)}s_{_N}(n) x(n)M×1​​=a(θ1​)s1​​(n)+a(θ2​)s2​​(n)+⋯+a(θN​)sN​​(n)利用小学二年级学过的矩阵知识:)，可以整理为 x ( n ) M × 1 = [ a ( θ 1 ) , a ( θ 2 ) , ⋯   , a ( θ N ) ] M × N × [ s 1 ( n ) , s 2 ( n ) , ⋯   , s N ( n ) ] 1 × N T ≜ A M × N s N × 1 \mathbf{x}(n) _{_{M\times 1}}= [\mathbf{a(\theta_1)},\mathbf{a(\theta_2)},\cdots,\mathbf{a(\theta_N)}]_{_{M\times N}}\times [s_{_1}(n), s_{_2}(n),\cdots, s_{_N}(n)]^T_{_{1\times N}}\triangleq\mathbf{A_{_{M\times N}}s_{_{N\times 1}}} x(n)M×1​​=[a(θ1​),a(θ2​),⋯,a(θN​)]M×N​​×[s1​​(n),s2​​(n),⋯,sN​​(n)]1×N​T​≜AM×N​​sN×1​​ A M × N \mathbf{A}_{M\times N} AM×N​也称为方向矩阵、导向矩阵等等。和导向矢量是一个意思，以后统称为导向矢量。

这里 A M × N \mathbf{A}_{M\times N} AM×N​的维度值得关注：

注意到接收信号的形式 x ( n ) = [ 1 , e − j 2 π f 0 d s i n ( θ ) c , e − j 2 π f 0 2 d s i n ( θ ) c , e − j 2 π f 0 3 d s i n ( θ ) c , ⋯   , e − j 2 π f 0 ( M − 1 ) d s i n ( θ ) c ] s ( n ) \mathbf{x}(n) = [1,e^{-j2\pi f_0\frac{dsin(\theta)}{c}},e^{-j2\pi f_0\frac{2dsin(\theta)}{c}},e^{-j2\pi f_0\frac{3dsin(\theta)}{c}},\cdots,e^{-j2\pi f_0\frac{(M-1)dsin(\theta)}{c}}]s(n) x(n)=[1,e−j2πf0​cdsin(θ)​,e−j2πf0​c2dsin(θ)​,e−j2πf0​c3dsin(θ)​,⋯,e−j2πf0​c(M−1)dsin(θ)​]s(n)虽然信号的角度我们不知道，但还是对于给定的阵列，其导向矢量的数学形式是知道的。比如对于ULA，肯定是Vandermonde结构的。就凭这一点，我们就有了一种DOA估计的方法。

具体而言，我们可以构造出一个导向矢量来，其中的角度任意给定（盲猜，不妨设为 α \alpha α）,那么就可以构造出一个来波方向为 α \alpha α的导向矢量为 a ( α ) = [ 1 , e − j 2 π f 0 d s i n ( α ) c , e − j 2 π f 0 2 d s i n ( α ) c , e − j 2 π f 0 3 d s i n ( α ) c , ⋯   , e − j 2 π f 0 ( M − 1 ) d s i n ( a l p h a ) c ] \mathbf{a(\alpha)}=[1,e^{-j2\pi f_0\frac{dsin(\alpha)}{c}},e^{-j2\pi f_0\frac{2dsin(\alpha)}{c}},e^{-j2\pi f_0\frac{3dsin(\alpha)}{c}},\cdots,e^{-j2\pi f_0\frac{(M-1)dsin(\\alpha)}{c}}] a(α)=[1,e−j2πf0​cdsin(α)​,e−j2πf0​c2dsin(α)​,e−j2πf0​c3dsin(α)​,⋯,e−j2πf0​c(M−1)dsin(alpha)​]用我们盲猜的这个导向矢量 a ( α ) \mathbf{a(\alpha)} a(α)和接收信号做向量内积，即 y = a H ( α ) ⋅ x ( n ) = a H ( α ) a ( θ ) s ( n ) y=\mathbf{a^H(\alpha)}\cdot\mathbf{x}(n)=\mathbf{a^H(\alpha)}\mathbf{a(\theta)}s(n) y=aH(α)⋅x(n)=aH(α)a(θ)s(n)结果应该是一个标量。简单计算一下可以得到 y = [ 1 + e j 2 π f 0 d s i n ( α ) − s i n ( θ ) c + e j 2 π f 0 d 2 [ s i n ( α ) − s i n ( θ ) ] c + ⋯ + e j 2 π f 0 d ( M − 1 ) [ s i n ( α ) − s i n ( θ ) ] c ] s ( n ) ≤ M s ( n ) y=[1+e^{j2\pi f_0d\frac{sin(\alpha)-sin(\theta)}{c}}+e^{j2\pi f_0d\frac{2[sin(\alpha)-sin(\theta)]}{c}}+\cdots+e^{j2\pi f_0d\frac{(M-1)[sin(\alpha)-sin(\theta)]}{c}}]s(n)\leq Ms(n) y=[1+ej2πf0​dcsin(α)−sin(θ)​+ej2πf0​dc2[sin(α)−sin(θ)]​+⋯+ej2πf0​dc(M−1)[sin(α)−sin(θ)]​]s(n)≤Ms(n)等号在 α = θ \alpha=\theta α=θ 处取得。

通过这个不等式可知，如果我们盲猜对了，即 α = θ \alpha=\theta α=θ，那么得到的结果是一个最大值。 因此，我们可以把所有的角度都猜一遍，找出那么结果结果最大的，其对应的角度就是我们DOA估计的结果。

这里就可以引出一个DOA估计的方法：

for: θ \theta θ from − 90 ° -90\degree −90° to − 90 ° -90\degree −90° caculate y α = a H ( α ) ⋅ x ( n ) y_{_\alpha}=\mathbf{a^H(\alpha)}\cdot\mathbf{x}(n) yα​​=aH(α)⋅x(n); end for θ = max ⁡ α y α \theta= \max\limits_{\alpha} y_{_\alpha} θ=αmax​yα​​

实例1：假设只有一个目标位于 θ 1 = 5 ° \theta_1=5\degree θ1​=5°的地方，得到的结果为 实例2：假设两个目标分别位于 θ 1 = 5 ° , θ 2 = 10 ° \theta_1=5\degree,\theta_2=10\degree θ1​=5°,θ2​=10°,采用上述方法得到的结果为实例3：假设两个目标分别位于 θ 1 = 5 ° , θ 2 = 30 ° \theta_1=5\degree,\theta_2=30\degree θ1​=5°,θ2​=30°,采用上述方法得到的结果为

从三个仿真实例可以看出，单目标没有问题，但是当两目标过于靠近时，此DOA算法不能分辨两个目标。这就带来了几个思考：

这些问题以后会慢慢说明。

代码已上传，有需要可点击此处下载

PS：原创，纯手敲。如有错误、建议等，请留言或私信，大家共同进步。